#include <iostream>
#include "includes/fv-evolve.hh"
#include "includes/slope.hh"
FVEvolve::FVEvolve(){
};
FVEvolve::~FVEvolve(){
};
void FVEvolve::fv_evolve(SubprojectBase* subproject, double* globalcoord, double* unitouternormal, double pillarlen, int i, int j, int currentdim, int* interscounter, double rdt){
  double unitouternormalsum=0.0;
  double var_i[memory::phynumber];
  double var_j[memory::phynumber];
  double fequ_i[memory::phynumber];
  double fequ_j[memory::phynumber];
  double slope_i[memory::phynumber];
  double slope_j[memory::phynumber];
  double var_left[memory::phynumber];
  double var_right[memory::phynumber];  
  double fequ_left[memory::phynumber];
  double fequ_right[memory::phynumber];  
  double flux[memory::phynumber];
  double dummyflag=0.0;
  for (int tmpi=0; tmpi<gridsetting::dimension; tmpi++)
    unitouternormalsum+=unitouternormal[tmpi];
  for (int tmpi=0; tmpi<memory::phynumber; tmpi++){
    var_i[tmpi]=memory::variable_tmp[tmpi][i];
    flux[tmpi]=0.0;
    var_j[tmpi]=0.0;
    fequ_i[tmpi]=0.0;
    fequ_j[tmpi]=0.0;    
    slope_i[tmpi]=0.0;
    slope_j[tmpi]=0.0;
    var_left[tmpi]=0.0;
    var_right[tmpi]=0.0;
    fequ_left[tmpi]=0.0;
    fequ_right[tmpi]=0.0;
  }
  // for first-order algorithem
  if (j>=0) {
    for(int tmpi=0; tmpi<memory::phynumber; tmpi++)
      var_j[tmpi]=memory::variable_tmp[tmpi][j];
  } else {
    subproject->boundary->boundary_function(globalcoord,i,unitouternormalsum,var_j,dummyflag);
  }
  subproject->equation->get_equation_flux(gridsetting::dimension,unitouternormal,var_i,fequ_i);
  subproject->equation->get_equation_flux(gridsetting::dimension,unitouternormal,var_j,fequ_j);
  // for second-order algorithem
  for (int tmpi=0; tmpi<memory::phynumber; tmpi++){
    slope_i[tmpi]=memory::limit_slope[tmpi][currentdim][i];
    slope_j[tmpi]=memory::limit_slope[tmpi][currentdim][j];  
    if (unitouternormal[currentdim]>0) {
      var_left[tmpi]=var_i[tmpi]+slope_i[tmpi]*pillarlen/2.0;
      var_right[tmpi]=var_j[tmpi]-slope_j[tmpi]*pillarlen/2.0;
    } else {
      var_left[tmpi]=var_j[tmpi]+slope_j[tmpi]*pillarlen/2.0;
      var_right[tmpi]=var_i[tmpi]-slope_i[tmpi]*pillarlen/2.0;
    }
  }
  subproject->equation->get_equation_flux(gridsetting::dimension,unitouternormal,var_left,fequ_left);
  subproject->equation->get_equation_flux(gridsetting::dimension,unitouternormal,var_right,fequ_right);  
  
  // --- simplest algorithm ------
  switch(setting::spacediscrete){
  case setting::lf:
    for(int tmpi=0; tmpi<memory::phynumber; tmpi++)
      flux[tmpi] = (fequ_i[tmpi]+fequ_j[tmpi])/2.0+unitouternormalsum*pillarlen/(rdt*2.0*gridsetting::dimension)*(var_i[tmpi]-var_j[tmpi]);
    break;
  case setting::upwind:
    // TODO: muscl
    for(int tmpi=0; tmpi<memory::phynumber; tmpi++)
      if (unitouternormalsum<0) {
	flux[tmpi] = std::max(fequ_j[tmpi],0.0)+std::min(fequ_i[tmpi],0.0);
      }else {
	flux[tmpi] = std::max(fequ_i[tmpi],0.0)+std::min(fequ_j[tmpi],0.0);
      }
    break;
  case setting::lw:
    // TODO: nonlinear, equ::u0*
    for(int tmpi=0; tmpi<memory::phynumber; tmpi++)
      flux[tmpi] = (fequ_i[tmpi]+fequ_j[tmpi])/2.0+unitouternormalsum/pillarlen*rdt/(2.0*gridsetting::dimension)*(fequ_i[tmpi]-fequ_j[tmpi]); 
    break;
  }
  // --- 2-order algorithm ------  
  switch(setting::spacediscrete){  
  case setting::tvdlf:
    // std::cout<<"fv-evolve.cc pillarlen: "<<pillarlen<<std::endl;        
    for(int tmpi=0; tmpi<memory::phynumber; tmpi++)
      flux[tmpi] = (fequ_i[tmpi]+fequ_j[tmpi])/2.0-0.5*pillarlen/rdt*timerule::cfl/gridsetting::dimension*(var_right[tmpi]-var_left[tmpi]);
    // flux[tmpi] = (fequ_i[tmpi]+fequ_j[tmpi])/2.0-0.5*(var_right[tmpi]-var_left[tmpi]);    
    //  flux[tmpi] = (fequ_left[tmpi]+fequ_right[tmpi])/2.0-0.5*rdt/pillarlen*(var_right[tmpi]-var_left[tmpi]);    
    // std::cout<<"fv-evolve.cc rdt/pillarlen: "<<rdt/pillarlen<<std::endl;            
    break;
  // case setting::richtmyerlw:
  //   for(int tmpi=0; tmpi<memory::phynumber; tmpi++){
  //     if(timerule::round_number==0) // first round is  LF
  //       flux[tmpi] = (fequ_i[tmpi]+fequ_j[tmpi])/2.0+unitouternormalsum*pillarlen/(rdt*2.0*gridsetting::dimension)*(var_i[tmpi]-var_j[tmpi]);
  //     if(timerule::round_number==1) // scecond round is center
  //   	flux[tmpi] = 0.0;
  //   }
  //   break;
  case setting::muscl:
    break;
  case setting::fct:
    break;
  }
  // --- weno algorithm ------

  // caculate update

  for(int tmpi=0; tmpi<memory::phynumber; tmpi++)
    memory::update[timerule::roundnumber][tmpi][i] -= 1.0/pillarlen*unitouternormalsum*flux[tmpi];  
};
